clear
clc

%%%%%%%%%%%%%%%
% Import Data % 
%%%%%%%%%%%%%%%

load('main_data_supply.mat');

for i_k = 1:4092

i_k

k = index_price_coeff(i_k,1);   


data_k = data_summary{k};
numProdsTotal = size(data_k,1);
X_MPEC_k = X_MPEC_summary{k};
X_MPEC_opt = X_MPEC_opt_summary{k};
share_k = data_k(:,5);
price_k = data_k(:,7) + 1;
gdp_capita_k = data_k(:,13);
distance_k = data_k(:,14);
population_k = data_k(:,16);
population_log_k = log(population_k);
Y_k = Y_summary{k};
FE_k = FE_summary{k};
market_id_k = data_k(:,end);
T = market_id_k(end,1);
nn = size(Y_k,2);
alpha_constant = X_MPEC_opt(1,1);
alpha_distance = X_MPEC_opt(2,1);
alpha_population = X_MPEC_opt(3,1);
alpha_price = X_MPEC_opt(4,1);
alpha_1 = X_MPEC_opt(12,1);
alpha_FE = X_MPEC_opt((3*(K_A+1) + 1):(3*(K_A+1) + (K - K_A)), 1);

xi_k = X_MPEC_opt((3*(K_A+1) + (K - K_A) + 2 + 1):(3*(K_A+1) + (K - K_A) + 2 + numProdsTotal), 1) + 10;     % Ensure positive values
xi_k = xi_k - alpha_constant - log(price_k)*alpha_price - alpha_population*population_log_k - alpha_distance*distance_k;
q_k = alpha_constant + xi_k + FE_k*alpha_FE;


prodsMarket = zeros(T,1);
prodsMarket_temp = data_k(:,17);
for t = 1:T
prodsMarket_temp2 = prodsMarket_temp(market_id_k==t);
prodsMarket(t,1) = prodsMarket_temp2(1,1);
end

marketStarts = zeros(T,1);
marketEnds = zeros(T,1);
marketStarts(1) = 1;
marketEnds(1) = prodsMarket(1);
for t=2:T
    marketStarts(t) = marketEnds(t-1) + 1;
    marketEnds(t) = marketStarts(t) + prodsMarket(t) - 1;
end

marketForProducts = zeros(numProdsTotal,1);
for t=1:T
marketForProducts(marketStarts(t):marketEnds(t)) = t;
end
numProdsTotal = size(data_k, 1);


delta_k = alpha_constant + xi_k + log(price_k)*alpha_price + alpha_population*population_log_k + alpha_distance*distance_k;
for t = 1:T   
Y_market = Y_k(t,:);
expmu(marketStarts(t):marketEnds(t),:) = exp(log(price_k(marketStarts(t):marketEnds(t),:))*alpha_1*Y_market + repmat(FE_k(marketStarts(t):marketEnds(t),:)*alpha_FE,1,nn));
mu_true(marketStarts(t):marketEnds(t),:) = log(price_k(marketStarts(t):marketEnds(t),:))*alpha_1*Y_market + repmat(FE_k(marketStarts(t):marketEnds(t),:)*alpha_FE,1,nn);
end
expmeanval = exp(delta_k);
meanval_true = delta_k;
oo = ones(1,nn);                  
sharesum = sparse(zeros(T,numProdsTotal));  % used to create denominators in logit predicted shares (i.e. sums numerators)
for t = 1:T
    sharesum(t,marketStarts(t):marketEnds(t)) = 1;
end
numer = (expmeanval*oo ).*expmu; 
v_true = (meanval_true*oo ) + mu_true; 

%==========================================================================


% Compute Utility and Price Index
%==========================================================================

% Factual
%--------------------------------------------------------------------------

filename_costs = sprintf('cost_estimates/costs_%i_k.mat', i_k);  
load(filename_costs, 'price_new', 'q_new', 'q_k', 'N_true', 'N_new', 'f_flag', 'f_fval', 'mc_opt_summary', 'p_flag', 'p_fval')
price_k_old = price_k;
price_k = price_new;
q_k_old = q_k;
q_k = q_new;

delta_k = alpha_constant + xi_k + log(price_k)*alpha_price + alpha_population*population_log_k + alpha_distance*distance_k;
delta_k = delta_k + (q_k - q_k_old) + log(N_new) - log(N_true);

for t = 1:T   
Y_market = Y_k(t,:);
mu_counter(marketStarts(t):marketEnds(t),:) = log(price_k(marketStarts(t):marketEnds(t),:))*alpha_1*Y_market + repmat(FE_k(marketStarts(t):marketEnds(t),:)*alpha_FE,1,nn);
expmu(marketStarts(t):marketEnds(t),:) = exp(log(price_k(marketStarts(t):marketEnds(t),:))*alpha_1*Y_market + repmat(FE_k(marketStarts(t):marketEnds(t),:)*alpha_FE,1,nn));
end
meanval_counter = delta_k;
expmeanval = exp(delta_k);
oo = ones(1,nn);                  
sharesum = sparse(zeros(1,size(data_k,1)));  % used to create denominators in logit predicted shares (i.e. sums numerators)
sharesum = sparse(zeros(T,numProdsTotal));  % used to create denominators in logit predicted shares (i.e. sums numerators)
for t = 1:T
    sharesum(t,marketStarts(t):marketEnds(t)) = 1;
end
numer_counter = (expmeanval*oo ).*expmu;  
v_counter = (meanval_counter*oo ) + mu_counter; 


P_true = zeros(T,nn);
beta_price = Y_k * alpha_1 + alpha_price;
for tt = 1:T
v_true_t = v_true(marketStarts(tt,1):marketEnds(tt,1),:,:); 
beta_price_t = beta_price(tt, :);
exp_v_true_t = exp(v_true_t);
s_exp_v_true_t = sum(exp_v_true_t);
ln_s_exp_v_true_t = log(s_exp_v_true_t);
P_true(tt,:) = exp(1) * exp(1./beta_price_t .* ln_s_exp_v_true_t) .* gamma(1 - 1./beta_price_t);
end

% Counterfactual
%--------------------------------------------------------------------------

P_counter = zeros(T,nn);
beta_price = Y_k * alpha_1 + alpha_price;
for tt = 1:T
data_k_t = data_k(marketStarts(tt,1):marketEnds(tt,1),:);   
index_domestic = data_k_t(:,3) == data_k_t(:,4);
v_counter_t = v_counter(marketStarts(tt,1):marketEnds(tt,1),:,:); 
v_counter_t = v_counter_t(index_domestic, :);
beta_price_t = beta_price(tt, :);
exp_v_counter_t = exp(v_counter_t);
s_exp_v_counter_t = exp_v_counter_t;
ln_s_exp_v_counter_t = log(s_exp_v_counter_t);
if sum(index_domestic) == 1
P_counter(tt,:) = exp(1) * exp(1./beta_price_t .* ln_s_exp_v_counter_t) .* gamma(1 - 1./beta_price_t);
end
if sum(index_domestic) == 0
P_counter(tt,:) = exp(1) * ones(1, nn);
end
end


% Sort results
%--------------------------------------------------------------------------

beta_price = Y_k * alpha_1 + alpha_price;
ratio = P_counter ./ P_true;

ratio_sorted = zeros(T, nn);
P_true_sorted = zeros(T, nn);
P_counter_sorted = zeros(T, nn);
beta_price_sorted = zeros(T, nn);
Y_sorted = zeros(T, nn);
for tt = 1:T
[Y_k_tt, ind] = sort(Y_k(tt,:));
Y_sorted(tt,:) = Y_k_tt;
P_true_sorted(tt,:) = P_true(tt,ind);
P_counter_sorted(tt,:) = P_counter(tt,ind);
ratio_sorted(tt,:) = ratio(tt,ind);
beta_price_sorted(tt,:) = beta_price(tt,ind);
end


diff = ratio_sorted(:,20) - ratio_sorted(:,80);
P_true_sorted_10 = P_true_sorted(:,10);
P_true_sorted_20 = P_true_sorted(:,20);
P_true_sorted_80 = P_true_sorted(:,80);
P_true_sorted_90 = P_true_sorted(:,90);
P_counter_sorted_10 = P_counter_sorted(:,10);
P_counter_sorted_20 = P_counter_sorted(:,20);
P_counter_sorted_80 = P_counter_sorted(:,80);
P_counter_sorted_90 = P_counter_sorted(:,90);
beta_price_20 = beta_price_sorted(:,20);
beta_price_80 = beta_price_sorted(:,80);

%==========================================================================


% Save Output to file
%==========================================================================

year_t = zeros(T,1);
quarter_t = zeros(T,1);
declarant_t = zeros(T,1);
price_domestic = zeros(T,1);
gdp_capita_domestic = zeros(T,1);
domestic_producer = zeros(T,1);
for tt = 1:T
data_tt = data_k(marketStarts(tt,1):marketEnds(tt,1),:);     
year_t(tt,1) = data_tt(1,1);
quarter_t(tt,1) = data_tt(1,2);
declarant_t(tt,1) = data_tt(1,3);
index_domestic = data_tt(:,3) == data_tt(:,4);
if sum(index_domestic) == 1
price_domestic(tt,1) = data_tt(index_domestic,7);
gdp_capita_domestic(tt,1) = data_tt(index_domestic,13);
domestic_producer(tt,1) = 1;
end
end

P_counter_sorted = P_counter_sorted .* (imag(P_counter_sorted) == 0);
P_counter_sorted(isnan(P_counter_sorted)) = 0;
P_counter_sorted(P_counter_sorted == Inf) = 0;


P_counter_sorted_10_real = P_counter_sorted_10 .* (imag(P_counter_sorted_10) == 0);
P_counter_sorted_10_real(isnan(P_counter_sorted_10)) = 0;
P_counter_sorted_10_real(P_counter_sorted_10_real == Inf) = 0;
P_counter_sorted_10_real(isnan(P_counter_sorted_10_real)) = 0;
P_counter_sorted_20_real = P_counter_sorted_20 .* (imag(P_counter_sorted_20) == 0);
P_counter_sorted_20_real(isnan(P_counter_sorted_20)) = 0;
P_counter_sorted_20_real(P_counter_sorted_20_real == Inf) = 0;
P_counter_sorted_20_real(isnan(P_counter_sorted_20_real)) = 0;
P_counter_sorted_80_real = P_counter_sorted_80 .* (imag(P_counter_sorted_80) == 0);
P_counter_sorted_80_real(isnan(P_counter_sorted_80)) = 0;
P_counter_sorted_80_real(P_counter_sorted_80_real == Inf) = 0;
P_counter_sorted_80_real(isnan(P_counter_sorted_80_real)) = 0;
P_counter_sorted_90_real = P_counter_sorted_90 .* (imag(P_counter_sorted_90) == 0);
P_counter_sorted_90_real(isnan(P_counter_sorted_90)) = 0;
P_counter_sorted_90_real(P_counter_sorted_90_real == Inf) = 0;
P_counter_sorted_90_real(isnan(P_counter_sorted_90_real)) = 0;
q_new_real = q_new .* (imag(q_new) == 0);
q_new_real(q_new_real == Inf) = 0;
q_new_real(isnan(q_new)) = 0;
price_new_real = price_new .* (imag(price_new) == 0);
p_fval_real = p_fval;
p_fval_real(:,1) = p_fval_real(:,1) .* (imag(p_fval_real(:,1)) == 0);
p_fval_real(:,2) = p_fval_real(:,2) .* (imag(p_fval_real(:,2)) == 0);

filename_table = sprintf('price_index/product_specific_gini/diff_full_adj_table_%i_k.csv', i_k); 
%data_save = array2table([data_k(:,1:4), q_k_old, q_new_real, price_k_old, price_new_real, repmat(f_flag, numProdsTotal, 1), [f_fval; repmat(9, numProdsTotal - size(f_fval,1), 1)], mc_opt_summary(:,3), mc_opt_summary(:,2), p_flag, p_fval_real, P_true_sorted_10(marketForProducts,:), P_true_sorted_20(marketForProducts,:), P_true_sorted_80(marketForProducts,:), P_true_sorted_90(marketForProducts,:), P_counter_sorted_10_real(marketForProducts,:), P_counter_sorted_20_real(marketForProducts,:), P_counter_sorted_80_real(marketForProducts,:), P_counter_sorted_90_real(marketForProducts,:), beta_price_20(marketForProducts,:), beta_price_80(marketForProducts,:), mc_opt_summary(:,1)]);
data_save = array2table([data_k(:,1:4), q_k_old, q_new_real, price_k_old, price_new_real, f_flag, [f_fval; repmat(9, numProdsTotal - size(f_fval,1), 1)], mc_opt_summary(:,3), mc_opt_summary(:,2), p_flag, p_fval_real, P_true_sorted_10(marketForProducts,:), P_true_sorted_20(marketForProducts,:), P_true_sorted_80(marketForProducts,:), P_true_sorted_90(marketForProducts,:), P_counter_sorted_10_real(marketForProducts,:), P_counter_sorted_20_real(marketForProducts,:), P_counter_sorted_80_real(marketForProducts,:), P_counter_sorted_90_real(marketForProducts,:), beta_price_20(marketForProducts,:), beta_price_80(marketForProducts,:), mc_opt_summary(:,1), Y_sorted(marketForProducts,:), P_true_sorted(marketForProducts,:), P_counter_sorted(marketForProducts,:), alpha_price*ones(numProdsTotal, 1), alpha_1*ones(numProdsTotal, 1)]);
data_save.Properties.VariableNames(1:26) = {'Year' 'Quarter' 'Declarant' 'Partner', 'q_old', 'q_new', 'p_old', 'p_new', 'f_flag', 'f_fval', 'mc_flag', 'mc_fval', 'p_flag', 'p_fval1', 'p_fval2', 'P_true_sorted_10', 'P_true_sorted_20', 'P_true_sorted_80', 'P_true_sorted_90', 'P_counter_sorted_10', 'P_counter_sorted_20', 'P_counter_sorted_80', 'P_counter_sorted_90', 'beta_price_20', 'beta_price_80', 'mc_opt'};  
data_save.Properties.VariableNames((end-1):end) = {'alpha_price', 'alpha_1'};
data_save = data_save(data_save.Declarant == data_save.Partner, :);
writetable(data_save, filename_table, 'Delimiter',',','QuoteStrings',true);

clearvars -except data_summary FE_summary income_mu_sigma index_price_coeff K K_A KK objective_MPEC_summary price_coeff_MPEC price_table status_MPEC_index status_MPEC_summary v_summary X_MPEC_opt_summary X_MPEC_summary Y_summary

end

